Self-organization with equilibration: 
a model for the intermediate phase in rigidity percolation 
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Recent experimental results for covalent glasses suggest the existence of an intermediate phase 
attributed to the self-organization of the glass network resulting from the tendency to minimize its 
internal stress. However, the exact nature of this experimentally measured phase remains unclear. 
We modify a previously proposed model of self-organization by generating a uniform sampling of 
stress-free networks. In our model, studied on a diluted triangular lattice, an unusual intermediate 
phase appears, in which both rigid and floppy networks have a chance to occur, a result also observed 
in a related model on a Bethe lattice by Barre et al. [Phys. Rev. Lett. 94, 208701 (2005)]. Our 
results for the bond-configurational entropy of self-organized networks, which turns out to be only 
about 2% lower than that of random networks, suggest that a self-organized intermediate phase 
could be common in systems near the rigidity percolation threshold. 

PACS numbers: 05.65.-|-b, 65.40.Gr, 61.43.Bn, 64.70.Pf 
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I. INTRODUCTION 

Rigidity theory 0,0, USB has improved considerably 
our understanding of the structural, elastic and dynam- 
icalproperties of systems such as chalcogenide glasses 
P, U 13 , interfaces Q and proteins , as a function of 
their connectivity. In its classical version, it introduces 
the concept of a rigidity transition, separating a soft (or 
floppy) and a rigid phases characterized by a mean coor- 
dination number. In many systems, optima or thresholds 
of various physical quantities are often observed at the 
rigidity transition. Some time ago, for example, Phillips 
[3 noted that among chalcogenides, the best glassform- 
ers have a mean coordination such that the number of 
degrees of freedom is equal to the number of covalent 
(bond-stretching and bond-bending) constraints. At this 
point, corresponding to the rigidity transition, networks 
are largely rigid but stress-free. This prevents crystalliza- 
tion for both kinetic and thermodynamic reasons: being 
stress-free, the glassy state is not too energetically un- 
favorable compared to the crystal; being rigid, the net- 
works lack the flexibility to efficiently explore the phase 
space and reach the crystalline state fast. Similarly, in 
the last 5 years, it has become clear that most proteins in 
their native state sit almost exactly at the rigidity transi- 
tion, which could be necessary to have enough flexibility 
to fuUflU their function while retaining their overall struc- 
ture 1^]. 

Recently, a series of experiments on chalcogenide and 
oxide glasses |lS|ll|lllillil|ll|llll3,IIllIl 
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Ha m m m m ig have demonstrated the exis- 
tence of a number of interesting and surprising behav- 
iors. For instance, glasses with near ly optim al prop- 
erties, such as the absence of aging El M m and 
vanishing non-reversing enthalpy of the glass transition 
liSIIllliMlillilliaSEllil are observed not just 
at a particular mean coordination, but in some range of 
coordinations, suggesting the presence of an intermediate 
phase between the floppy and the rigid phases. While the 
details and exact origin of this intermediate phase are still 
a matter of debate (see, for example, the recent experi- 
mental paper by Golovchak et al. 26]), it appears that it 
is due to the self-organization of the network minimizing 
the internal stress. If, according to Phillips' argument, 
we expect "optimal" glasses to be rigid but stress-free, 
we should now expect to find this property everywhere 
in the intermediate phase rather than only at a single 
critical point, as in the standard phase diagram of rigid- 
i ty^p erc olat ion ; there is now some direct evidence for this 

A few models were proposed to explain this self- 
organization. Thorpe et al. have shown in an 
out-of-equilibrium model that it is possible to generate 
a stress-free intermediate phase. Barre et al. have 
shown in addition that such a phase was thermodynam- 
ically stable on a ring-free Bethe lattice, where each site 
has 3 degrees of freedom, but the whole network is em- 
bedded in an infinite-dimensional space. Taking a dif- 
ferent approach, Micoulaut ^3 demonstrated that one 
could recover an intermediate phase by concentrating all 
the strain in local structures. 

The goal of this paper is to assess whether or not a 
self-organized network with a finite-dimensional topology 
is also thermodynamically possible. This verification is 
important on two counts: (1) the Bethe lattice is a loop- 
less structure producing a first-order rigidity transition 
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0, Isil Is^ while 2- and 3-dimensional regular networks 
undergo a second-order phase transition [J 1^; (2) the 
original model of Thorpe et al. is an out-of-equilibrium 
model which could lead to highly atypical networks. 

For simplicity, we study two-dimensional central-force 
(2D CF) networks. This allows us to consider networks 
of much bigger linear dimensions than is possible in 3D, 
limiting the finite-size effects. Moreover, in all cases stud- 
ied until now, rigidity results for 2D CF networks have 
been qualitatively very similar to those obtained for 3D 
glass networks. Our results should therefore also apply 
to glass networks. 

In the next section, we review the basic facts about 
rigidity and the intermediate phase. We then explain 
the model used here and present our results for a bond- 
diluted two-dimensional triangular network. We first dis- 
cuss an unusual property of our model: a network in the 
intermediate phase can be either rigid or floppy with a 
finite probability. We then focus on the calculation of 
the entropy of self-organized networks. 

II. THE INTERMEDIATE PHASE IN THE 
RIGIDITY PHASE DIAGRAM 

In rigidity theory, an elastic network is characterized 
by the number of motions, called floppy modes, that do 
not distort any constraints. In a system with no con- 
straints, all degrees of freedom are floppy modes and thus 
their number is dN , where N is the number of atoms and 
d is the dimensionality of space. In an approximation due 
to Maxwell and known as Maxwell counting ^| , it is as- 
sumed that each additional constraint removes a floppy 
mode, so that the number of floppy modes for a given 
number of constraints can be written as 

F = dN - Nc. (1) 

When the number of constraints becomes equal to the 
number of degrees of freedom, F — and the network 
undergoes a rigidity percolation transition, going from 
floppy to rigid. In a 2D CF network, the number of 
constraints per atom is (r)/2, where (r) is the mean 
coordination (the average number of connections of a 
site); the critical coordination is therefore (r) = 4. In 
chalcogenide glasses, characterized by the chemical for- 
mula Ax^yCi-x-y, where A is an atom of valence 4 (usu- 
ally Ge or Si), B is an atom of valence 3 (As or P) and 
C is a chalcogen (Se, S or Te), counting both covalent 
bonds and their related angular constraints, the critical 
coordination is (r) = 2.4 

Since Maxwell counting is a mean-field theory, it ig- 
nores fluctuations and correlations that can be built in 
the network. An overall rigid network can have some in- 
ternal localized floppy modes. Also, a constraint inserted 
into a piece of the network that is already rigid does not 
remove floppy modes. Such constraints, known as re- 
dundant, are obviously present in the rigid phase, where 
Eq. ^ gives a negative number of floppy modes, since 



this number cannot be lower than d{d+ l)/2 (6 in 3D, 3 
in 2D), to account for rigid body translations and rota- 
tions; but redundant constraints can also be present, of 
course, in an overall floppy network. In generic networks, 
such as glasses, where the bond lengths vary, these con- 
straints create an internal stress and increase the elastic 
energy, since at least some part of the network has to be 
deformed to accommodate them. 

If the number of redundant constraints, Nji, is known, 
then the Maxwell counting formula can be corrected: 

F^dN-N, + Nr. (2) 

This result is exact — the problem lies in calculating Nr . 
A theorem by Laman [s^ makes this possible. Consider 
all possible subnetworks of a system. If the number of de- 
grees of freedom minus the number of constraints is less 
than d{d + l)/2 for at least one subnetwork, there must 
be redundant constraints. Laman showed that this is the 
only way redundant constraints can occur: for them to 
be present, the above must be satisfied for at least one 
subnetwork. This is strictly true in 2D; although there 
are known counterexamples for general networks in 3D, 
it is assumed true (there is no rigorous mathematical 
proof, but no known counterexamples either) for glassy 
networks with covalent bonding including angular con- 
straints mni. 

Laman's theorem is used in a computer algorithm for 
rigidity analysis, known as the pebble game ,3, .37, .38,. .3^ 
1401 . The pebble game starts with an empty network and 
then one constraint is added at a time and each added 
constraint is checked for redundancy. Thus at every stage 
in the network-building process the number of redundant 
constraints, Nr, and, according to Eq. |21 F, are known 
exactly. Independent constraints are matched to pebbles 
that are assigned to sites and whose total number is equal 
to the number of degrees of freedom; thus the number of 
free pebbles (not matched to any constraint) gives the 
number of floppy modes. The pebble game, in addition, 
identifies stressed regions in the network and also offers 
rigid cluster decomposition that identifies all rigid clus- 
ters in the network. The analysis provided by the pebble 
game is purely topological, the details of the geometry 
are not taken into account, nor is the exact expression 
for the forces. The downside is that the pebble game 
(as well as Laman's theorem itself) is only applicable to 
generic networks: networks special in some way (having 
parallel bonds, for instance) may have rigidity proper- 
ties that are different from those of the vast majority of 
networks of a given topology. The fraction of such spe- 
cial (or non-generic) networks among all networks of a 
given topology is, however, zero, and so a covalent glass 
network, being disordered, can be safely assumed to be 
generic. 

Rece ntly , a series of experiments [Tol Hit IT^ IT^ Q, 
EllllllllllllSillllllllll have suggested that 
there could be not one but two phase transitions near 
(r) = 2.4, with the opening of an intermediate phase be- 
tween the phases already known. The properties of this 
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phase suggest that the intermediate phase is rigid but 
stress-free. To explain the presence of two transitions and 
the intermediate phase between them within the frame- 
work of rigidity theory, it has been proposed that 
the glass networks self- organize in some way. Broadly 
speaking, any reduction in the amount of disorder in the 
network as it tries to minimize its free energy can be re- 
ferred to as self-organization; chemical order, especially 
strong in oxide glasses like silica, would be one example. 
Thorpe et al. l23 | considered a particular kind of 
self-organization: glass networks minimizing their elastic 
stress energy. As a proof of principle, they constructed 
the following model. Starting with a low-coordination 
stress-free network, bonds are added one at a time with 
the restriction that they cannot be redundant and thus 
add stress to the network, using the pebble game for con- 
structing and analyzing the network at each step. This 
process is repeated until it is no longer possible to add 
a bond without introducing stress to the network. After 
that, bond insertion continues at random. The maximum 
coordination at which no stress is present, according to 
Eq. cannot exceed the rigidity threshold according 
to Maxwell counting ((r) =4 for 2D CF networks and 
(r) = 2.4 for covalent glass networks). In fact, in the 
model of Thorpe et al., the Maxwell counting threshold 
value is reached without stress for 2D CF networks, but 
not for the 3D glass network In this model, once the 
stress appears, it immediately percolates, corresponding 
to the upper boundary of the intermediate phase. In gen- 
eral, this need not be the case, since a network can have 
finite stressed regions without stress percolating (this is 
the case in the model of Micoulaut and Phillips jSfl EH)- 

To observe the rigidity transition in a two-dimensional 
diluted regular lattice, one needs a lattice with coordina- 
tion bigger than 4, and so the triangular lattice is a nat- 
ural choice. Since the pebble game algorithm can only 
be applied to generic networks, one has to assume that 
the triangular lattice is distorted (for example, by hav- 
ing some disorder in bond lengths). Previous numerical 
studies for the randomly diluted triangular lattice with- 
out self-organization indicate a single rigidity and stress 
transition at (r) = 3.961 — very close to the Maxwell 
counting prediction (Fig. ^ Q . Note that the rigidity 
and stress transitions coincide in this case. 

In the self-organization model of Thorpe et al., a nu- 
merical simulation reveals instead two phase transitions 
(Fig. ^ . As the coordination is increased from the floppy 
phase, a percolating rigid cluster appears at (r) — 3.905. 
At this point, by construction, the network is still stress- 
free. This is the lower boundary of the rigid but stress- 
free intermediate phase — the rigidity percolation transi- 
tion. As the mean coordination continues to grow, keep- 
ing the network stress-free becomes impossible. At this 
point, which, as mentioned above, is at (r) = 4, the stress 
appears and immediately percolates. This is the stress 
percolation transition, which is the upper boundary of 
the intermediate phase. 

The probability of a percolating cluster in the model 
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FIG. 1: The dependence of the fraction of bonds in the per- 
colating rigid cluster and the percolating stressed region for 
the random case and the self-organized case without equili- 
bration. In the random case, the rigidity and stress perco- 
lation thresholds coincide. In the self-organized case, these 
thresholds do not coincide and there is an intermediate phase 
(shaded) between them. The stressed region is defined here 
as the contiguous percolating set of stressed bonds within the 
percolating rigid cluster. The simulations are averaged over 2 
realizations on the bond-diluted triangular lattice of 400 x 400 
sites. The figure is taken from Ref. |2q| . 



of Thorpe et al. is zero below the rigidity percolation 
threshold and one above in the thermodynamic limit, as 
normally is the case for percolation transitions. This is 
illustrated in Fig. |21 where the probability of having a 
percolating cluster is shown for different network sizes. 
As expected, the dependence gets closer to a step func- 
tion as the size increases. 



III. A SELF-ORGANIZATION MODEL WITH 
EQUILIBRATION 

The self-organization model of Thorpe et al. is peculiar 
in that bonds are only added to the network and never 
removed. While one can imagine a very rapid quench pro- 
cess in which indeed bond formation dominates, this pro- 
cess does not lead to formation of good glasses. There- 
fore a way of building equilibrated stress-free networks is 
needed. As the elastic energy of a stress-free network is 
zero, any such networks should occur with equal proba- 
bility. 

A similar issue arose before in the case of conven- 
tional (or connectivity) percolation. In the rigidity case, 
self-organization proceeds by avoiding stress or redun- 
dancy, i.e., bonds that connect already mutually rigid 
sites. The connectivity analog consists in avoiding con- 
nections between sites that are already connected, i.e., 
creating loops. Straley ^3 proposed a model directly 
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FIG. 2: The fraction of networks in which the percolating 
rigid cluster is present as a function of (r) for the model of self- 
organization without equilibration. Each curve is obtained 
from 100 separate runs on a bond-diluted triangular lattice; 
the lattice sizes are indicated in the legend. 



analogous to the one by Thorpe et al, i.e., with bonds 
inserted one at a time and those forming loops rejected; 
connectivity percolation occurs at some point, and then 
there is an "intermediate phase" (although it was not re- 
ferred to as such in the connectivity percolation context) 
that is connected but without any loops, until a point is 
reached at which avoiding loops is no longer possible (at 
this point the network is a spanning tree). It was real- 
ized later on that this model does not produce a uniform 
ensemble, in which every loopless network would occur 
with equal probability Several authors, using a variety 
of methods 0i Ell j claimed then that in the equi- 
librated uniform ensemble, connectivity percolation does 
not occur until the spanning tree limit is reached, i.e., 
there is no intermediate phase. This shows that the re- 
sults may change significantly depending on the ensemble 
of self-organized networks that is considered. 

Braswell et al. 01 1 in particular, used the following al- 
gorithm to generate equiprobable loopless networks: take 
an arbitrary loopless network, choose a bond at random, 
delete it and then reinsert at an arbitrary place where it 
would not form a loop, with this place chosen equiproba- 
bly among all such places. They showed that after equi- 
libration has taken place, this method would indeed gen- 
erate the uniform ensemble of networks, by proving the 
detailed balance condition, i.e., that given some loopless 
network 1, the probability that in a single step of the al- 
gorithm some other network 2 would be produced is the 
same as the probability of going in the opposite direction, 
i.e., from network 2 to network 1. Their arguments fully 
apply to the rigidity case as well. 

In view of the above, we consider a variety of the self- 



organization model by Thorpe et at, adding equilibration 
that produces equiprobable stress-free networks. Like in 
the previous model, we start with the "empty" network 
without bonds and start inserting bonds one by one with- 
out creating stress. After every bond insertion, we equi- 
librate by doing bond swaps following the procedure de- 
scribed above, i.e., choose a bond at random, delete it and 
then insert a bond elsewhere choosing at random among 
the places where that new bond would not create stress. 
It is worth noting that in general it is rather difficult to 
handle removal of constraints within the pebble game; 
but it is easy to remove an unstressed constraint, as this 
simply involves releasing the associated pebble with no 
additional pebble rearrangement. Since in our case all 
constraints are unstressed by construction, no problem 
arises. In this paper we focus on the intermediate phase 
and do not investigate the stressed phase, so we stop at 
the point at which further stressless insertion becomes 
impossible. 

A very similar model has been proposed recently by 
Barre et al. using a Bethe lattice and, as an added 
sophistication, an energy-cost function linear in the num- 
ber of redundant constraints. While this energy is some- 
what unrealistic, it provides a thermodynamic justifica- 
tion for the existence of the intermediate phase. How- 
ever, Bethe lattices are particular constructions, leading 
to a first-order rigidity phase transition while 2D central- 
force and 3D bond-bending networks show a second-order 
rigidity transition. By comparison, our model in essence 
assigns an infinite cost to redundant bonds and corre- 
sponds therefore to the T ^ limit of Barre's model but 
on a regular lattice with a second-order rigidity transi- 
tion. 



IV. RIGIDITY PERCOLATION IN THE 
INTERMEDIATE PHASE 

Our simulations for the new model with equilibration 
are done for 2D CF bond-diluted triangular lattices. Pe- 
riodic boundary conditions were used with the supercell 
consisting of the same number of unit cells in both di- 
rections. We have chosen the duration of equilibration 
equal to 100 steps above (r) = 3.5 and 10 steps below 
(where even in random networks there are very few re- 
dundant constraints, so the self-organized networks are 
almost completely random anyway and long equilibra- 
tion is not needed). This is sufficient for convergence as 
shown in Fig. O for a 100 x 100 lattice; we also checked 
that 100 steps was sufficient by comparing with a very 
long equilibration run at a single point (r) =3.95 for the 
largest lattice used here (not shown). The result wc get 
is quite different from that obtained without equilibra- 
tion. Fig. 21 just as Fig. [3 shows the probability of rigid- 
ity percolation as a function of (r) for several different 
sizes, but now for the model with equilibration. The self- 
organization still opens an intermediate phase around the 
critical point found at (r) « 3.961 in the standard rigid- 
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FIG. 3: Same as in Fig. m but now for the model with equili- 
bration, for triangular networks of 10000 sites and for several 
different equilibration times indicated in the legend as the 
number of equilibration steps per each inserted bond above 
(r) — 3.5; below (r) = 3.5, there are 10 equilibration steps 
per bond in all cases. The data are likewise from 100 sep- 
arate runs, but in addition, from each run all networks ob- 
tained during the equilibration procedure at the given mean 
coordination are taken into account. 



ity phase diagram. But rather than approaching a step 
function as the size increases, the result for the proba- 
bility of percolation is now a gradual increase from at 
(r) « 3.94 to 1 at (r) = 4. This is similar to the result on 
a Bethe lattice in the model of Barre et ai, even though 
the rigidity percolation transition is of a different order. 
The dependence of the percolation probability on (r) is 
close to linear and this linearity may, in fact, be exact, 
although a very small non-linear region near the lower 
boundary of the intermediate phase cannot be ruled out. 
One difference between our result and that presented by 
Barre et al. is that since their consideration is at a non- 
zero temperature, there is always a possibility of having 
a small number of redundant constraints; since adding 
very few (perhaps C(l)) redundant constraints to an un- 
stressed network is often enough for stress percolation, it 
is not surprising that they have found a finite probability 
of both rigidity and stress percolation in the intermediate 
phase, whereas in our case, the stress percolation proba- 
bility is, of course, zero by construction. 

V. ENTROPY COST OF SELF-ORGANIZATION 

Although we do not include a potential energy explic- 
itly, the self-organization models discussed here are con- 
structed to prevent the build-up of stress, implicitly min- 
imizing the potential energy of the system. At a finite 
temperature, however, we are interested in minimizing 



FIG. 4: Same as in Fig. |21 but for different network sizes 
indicated in the legend, always using 10 equilibration steps 
per inserted bond below (r) — 3.5 and 100 equilibration steps 
above. 



the free energy. As glasses are formed at a non-zero 
temperature, the thermodynamical state will be influ- 
enced by the balance between the entropic cost associ- 
ated with generating a self-organized network and the 
energetic cost of creating internal stress: if the entropy of 
the random network is large compared with that of the 
self-organized one, then it is likely that very little self- 
organization will take place and the problem discussed 
here becomes irrelevant. 

The entropy of a covalent network can be viewed as 
consisting of two parts. The first part, the topological en- 
tropy, is proportional to the logarithm of the number of 
different possible bond topologies or bond configurations. 
The second part, which can be called, somewhat simplis- 
tically, the flexibility entropy, depends on the phase space 
available to each such bond configuration. This division 
of the total entropy is similar, but not identical, to the 
traditional division into the configurational and vibra- 
tional entropy in the inherent structure formalism |46j| . 
In particular, flexible networks exhibit a wide range of 
motions and would generally correspond to more than 
one inherent structure, when a potential energy function 
is defined. For this reason, the flexibility entropy includes 
both harmonic and anharmonic contributions associated 
with a given topology of the covalent network; its exact 
evaluation is difhcult and goes beyond the scope of this 
paper. However, since the flexibility entropy is expected 
to be roughly proportional to the number of floppy modes 
[irl lisL |49| , and, according to Eq. Q , the number of 
floppy modes in a self-organized network (with Nn = 0) 
is smaller than in a random network {Nn > 0) with the 
same number of constraints Nc, the flexibility entropy 
of the self-organized network is likely to be smaller than 
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that of the random network. This difference is proba- 
bly not very large, especially in real systems where long- 
range forces reduce significantly the available configura- 
tion space even in the floppy phase. 

The topological entropy is, in general, difficult to cal- 
culate as well, although methods, such as that by Vink 
and Barkema *5ff|, exist. It is much simpler for a lattice- 
based model like ours, as it requires only counting the 
number of possible bond configurations on a lattice. In 
this case, the topological entropy (which we can also call 
the bond-configurational entropy) is simply 



5bc((r)) =ln7Vbc((r)), 



(3) 



where Nbciif)) is the number of stress- free configurations 
with mean coordination (r) and the Boltzmann constant 
fcs is put equal to 1. To calculate A'bc, we use the fol- 
lowing approach. Suppose the number of stress-free net- 
works having Nb bonds, Ni^dNs), is known. From a 
stress-free network having Nb bonds, it is possible to 
produce a stress-free network having Nb + 1 bonds by 
adding a bond in one of those places where this added 
bond would not create stress. Suppose on average there 
are n+{NB) such places or ways to create a stress- free 
network with Nb + 1 bonds. On the other hand, for any 
stress-free network with Nb + 1 bonds, there are always 
exactly n-{NB + 1) = Nb + 1 ways of creating a stress- 
free network with Nb bonds by removing any one of the 
Nb + 1 bonds. Moreover, if a network with Nb + 1 bonds 
can be created from a network with Nb bonds by adding 
a bond, then the latter network can always be obtained 
from the former by removing that same bond and vice 
versa, so the process is reversible. Then the number of 
stress-free networks with Nb + 1 bonds is 



Ni,,{NB + l) = Afbc(iVB) 
= Nbc{NB) 



n+jNB) 
n^{NB + l) 
n+jNe) 
Nb + 1' 



(4) 



If n-^-{NB) is known for all Nb, then, using iVbc(O) = 1 as 
the initial condition, Eq. Q can be iterated to yield all 
^bc(^s)- In practice, n+(7Vs) are obtained numerically, 
by a sort of Monte Carlo procedure, where some of the 
places in the network where a bond is missing (compared 
to the full undiluted lattice) are tried and it is found in 
what fraction of such places addition of a bond would 
not create stress. In our simulations, we use 100 such at- 
tempts per network. The result is then averaged over the 
networks with a given number of bonds obtained during 
the equilibration procedure. 

Note that the same procedure can be repeated for 
the random case (without any self-organization). In this 
case, rt-|_(iVs) should count all ways to create a new net- 
work of Nb + 1 bonds (no matter stress-free or not) out 
of a network of Nb bonds, and there are as many ways 
to do that as there are places where a bond is missing 
(compared to the full lattice); thus for a random net- 
work, rfj^i^NB) = Nb/ — Nb, where Nb/ is the number 



of bonds in the full lattice. On the other hand, the analog 
of quantity n^{NB + 1), which we denote rf_{NB -I- 1), is 
still Nb + 1. Therefore, for the random network. 



K,(7Vb + 1) = K^Nb) 
= Kc{Nb) 



n;(A^B) 

nL{NB + l) 
Nbj-Nb 
Nb + 1 



(5) 



This can, of course, be used to obtain the analytical result 
for any A''^, which is simply a binomial coefficient. If we 
are interested in the difference AS between the entropies 
of the random and self-organized networks, this is simply 



ASINe 



In 



^b^c(^^B) 

N^c{Nb)' 



(6) 



The change of this difference when a bond is added is 

AS{Nb + 1) - AS{Nb) = 

'^bc(^B + 1) N^Nb) 



= In 
= In 



NI^Nb) A^bc(^B + 1) 
n'+iNB) Nb + 1 



Nb + 1 n+{NB) 
= -IniyiNB), 



(7) 



where v{Nb) = n^(NB)/n^{NB) is the average /raction 
of bonds whose insertion would not create stress among 
all missing bonds in the network. Note that v{Nb) is 
what is calculated directly by the Monte Carlo procedure 
described above. 

The entropy S is, of course, an extensive quantity, i.e., 
it is proportional to the network size (for big sizes). We 
can introduce the entropy per bond of the full lattice, 
s = S/NBf- In the thermodynamic limit, 

d(As) _ AS{Nb + 1)-AS{Nb) 



d{r) 



N, 



Bf 



(r) 



Nb + 1 

In u 



(r) 



Nb 



Nbj 



Nb + 1 



(r) 



Nb 



(8) 



where (r) 



is the mean coordination of a network with 



Nb 



Nb bonds, and since (r) = 2Nb/N, 

Nb 



d(As) 



N 



2N, 



■ In V. 



Bf 



(9) 



In the full triangular lattice, the number of bonds Nbj 
is three times the number of sites N , so we get 



d(As) 
d(r) 



■ In u. 



(10) 



The quantity v obtained numerically is plotted in 
Fig. 13 From this plot, it seems to go to zero as (5 = 



Ql , 1 , 1 , 1 , 1 , 1 

3 3.2 3.4 3.6 3.8 4 
Mean coordination (r) 

FIG. 5: The fraction v of "allowed bonds" (places where a 
bond can be inserted without creating stress) as a function of 
the mean coordination (r), for a triangular network of 50176 
sites, with equilibration time of 1000 steps per inserted bond 
above (r) = 3.5 and 10 steps below. The result is obtained 
by a Monte Carlo procedure described in the text; at each 
mean coordination, it is averaged over the networks obtained 
during equilibration. 



4 — (r) — > 0+ and appears to change linearly as a func- 
tion of b in this limit. In general, if in this limit v ~ 5™, 
then, according to Eq. H10() . 

s((5) = —(5 In (5 + regular part. (11) 

Figure El shows the entropy difference As calculated 
by iterating Eq. Q, with v obtained numerically after 
every bond addition. Since the simulations are done for 
networks of a finite size and with finite equilibration time, 
an extrapolation to the infinite size was done by fitting 
the simulation results to the function 

M(r);iV,r) = ^ + ^ + a((.)), (12) 
N T 

where r is the equilibration time (in equilibration steps 
per added bond) used above (r) = 3.5 (below 3.5, we al- 
ways use just 10 steps, for reasons explained above). As 
in different runs for different sizes the data were taken at 
slightly different points, linear interpolation was some- 
times done to obtain the values at the same (r) in all 
cases. In total, data for 182 [N, r) combinations were 
used, with N between 5476 and 50176 sites and r be- 
tween 20 and 1000 steps (up to 10000 steps for 10000 
sites). Data for higher N and r were assigned higher 
weights in the fit, since increasing TV and r decreases the 
amount of noise in the data (the latter because the results 
for V are averaged over all equilibration steps, and thus 
the more equilibration steps there are the smaller the er- 
ror in u). Function C{{r)) representing the asymptotic 
value of the entropy difference is also shown in Fig. 



Ql I I 1 I I I I I 1 I 

3.9 3.92 3.94 3.96 3.98 4 
Mean coordination <r) 

FIG. 6: The difference in the bond-configurational entropy 
between random and self-organized triangular networks, As, 
as a function of the mean coordination (r) , shown for several 
network sizes and equilibration times, as well as the asymp- 
totic value (i.e., the extrapolation to the infinite size and re- 
laxation time, as described in the text). Quantity v used to 
calculate As is obtained by the Monte Carlo procedure de- 
scribed in the text and then averaged over networks obtained 
during equilibration at the given (r). In the main plot, all 
thin lines are results for 50176 sites; the equilibration times 
are, top to bottom: 10 steps, 30 steps, 100 steps and 300 steps 
per added bond. The thick line is the asymptotic. In the in- 
set, all thin lines are results for 1000 equilibration steps per 
added bond. The sizes are, top to bottom: 5476 sites, 15376 
sites and 50176 sites. The thick line is again the asymptotic. 
The asymptotics are themselves extrapolated to (r) = 4, as 
described in the text. All equilibration times listed are for 
(r) > 3.5; for (r) < 3.5, 10 equilibration steps per added 
bond are always used. 



For a finite-size network, it is only possible to reach a 
point 3 bonds short of (r) = 4 without creating stress. 
For this reason, our simulations were stopped somewhere 
around (r) = 3.999. Taking into account Eq. Hll|l . we 
have fitted the asymptotic entropy difference C{6) be- 
tween (r) = 3.97 and 3.999 using the following function: 

C{6) ^ ao + aiSlnS + + asS^. (13) 

The fit is essentially perfect, and the obtained value of 
fli = 0.1636 is consistent with the value of 1/6 expected 
for m = 1, according to Eq. Hll|l . The fit is used to 
complete the curve in Fig. up to (r) = 4. The value of 
the entropy difference at (r) = 4 is oq = 0.0146. This is 
the biggest value of the entropy difference, but it is still 
small, only about 2% of the bond-configurational entropy 
of the random network (which at (r) — 4, when 2/3 of the 
bonds are present, is -[(2/3) ln(2/3) + (1/3) ln(l/3)] = 
0.6365...). 

In the above calculation of the entropy we explicitly 
use the fact that all stress-free networks with a given 
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number of bonds are equiprobable in our new model. A 
similar consideration for the old self-organization model 
without equilibration would be much more difficult. 

We note, finally, that according to Eq. there is 

a non-analyticity in the behavior of the entropy when 
(r) 4^, i.e., at the stress transition. We should also 
expect some very weak non-analyticity (perhaps a break 
or a cusp in a higher derivative) at the lower bound- 
ary of the intermediate phase (the rigidity transition) at 
(r) « 3.94, but it is too weak to be seen in our simulation 
results. 



VI. CONCLUSION 

We have considered a model of self-organization in elas- 
tic networks, adding an equilibration feature to the model 
previously considered by Thorpe et al. [23, HI] . In our 
model, we find an intermediate phase in the rigidity phase 
diagram, in which the fraction of networks in which rigid- 
ity percolates is between and 1 in the range of mean 
coordination between 3.94 and 4.0 for the bond-diluted 
triangular lattice, a resultqualitatively similar to that 
obtained by Barre et al. [23 in a closely related model 
on a Bethe lattice. 

Calculating the bond-configurational entropy of these 
self-organized networks, we find that it is only about 2% 
smaller than that of randomly-connected networks. Pro- 
vided that the flexibility entropy, which should reduce 
the stability of the intermediate phase, is not so sensitive 
to the self-organization, the intermediate phase is likely 
to be present in most systems with the right range of 
mean coordination. 

Our results support the current explanation of the 
intermediate phase in chalcogenide glasses. Self- 
organization might also be important in the dynamics 
of proteins as they have a coordination near the critical 
value. In a real material, it is likely that the intermediate 
phase is not perfectly stress-free. A most likely structure 
will therefore be mostly unstressed with overconstrained 
local regions, in a mixture of the model presented here 



and that introduced recently by Micoulaut et al. |30ll4ll |. 

The fact that there is now a possibility of having non- 
rigid networks in the intermediate phase does not inval- 
idate the concept of this phase as lacking both excessive 
flexibility and stress. Indeed, even though the network 
may technically be floppy, say, because, of a single floppy 
mode or a very small number of such modes spanning the 
whole network, for any practical purposes, there would 
be no difference between such a network and a fully rigid 
and unstressed one, especially when the existence of the 
neglected weaker interactions is taken into account. 

Knowing that self-organization can exist from a ther- 
modynamical point of view, there is still considerable 
work to do in order to fully understand the intermediate 
phase. Among the obvious future directions of this work 
we can mention: repeating the simulations described here 
for 3D bond-bending glass networks; getting a better idea 
of the geometry of self-organized networks, in particular, 
possible long-range correlations in them; evaluating the 
flexibility entropy effects using a particular potential en- 
ergy function. 

Finally, there is a suspicious discrepancy between our 
results reported here and those obtained for a similar 
model in connectivity percolation jl^, IS l45l |. Even 
though it is in principle possible that there is an interme- 
diate phase in the rigidity case but not in the connectivity 
case, this seems very unlikely. Perhaps the time has come 
to re-evaluate these old results — we certainly have the 
benefit of the much increased computational power. 
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